classdef Helper< handle
    properties (Constant)
        leaveOne = false %%要修改CalcAllErrors
    end
    methods(Static = true)
        % % Step1：把原始数据导出为mat文件
        function loadAllData(datapath, bands, isRefl, fileMatPath)
            %     fileMatPath  = 'F:\12份土壤\matlab\data';
            
            fileName = fullfile(datapath, '1-6_up.txt')
            roie = MyRoi.LoadBatchRoies(fileName, bands, isRefl);
            
            roiRange = 1:6;
            suffix = '_1';
            MyRoi.SaveAllRoi(roie, fileMatPath, roiRange, suffix)
            disp('done1')
            
            %     fileName = 'F:\12份土壤\hou\1-6_down.txt'
            fileName = fullfile(datapath, '1-6_down.txt')
            roie = MyRoi.LoadBatchRoies(fileName, bands, isRefl);
            roiRange = 1:6;
            suffix = '_2';
            MyRoi.SaveAllRoi(roie, fileMatPath, roiRange, suffix)
            disp('done2')
            
            %     fileName = 'F:\12份土壤\hou\7-12_up.txt'
            fileName = fullfile(datapath, '7-12_up.txt')
            roie = MyRoi.LoadBatchRoies(fileName, bands, isRefl);
            
            roiRange = 7:12;
            suffix = '_1';
            MyRoi.SaveAllRoi(roie, fileMatPath, roiRange, suffix)
            disp('done3')
            
            %     fileName = 'F:\12份土壤\hou\7-12_down.txt'
            fileName = fullfile(datapath, '7-12_down.txt')
            roie = MyRoi.LoadBatchRoies(fileName, bands, isRefl);
            
            roiRange = 7:12;
            suffix = '_2';
            MyRoi.SaveAllRoi(roie, fileMatPath, roiRange, suffix)
            disp('done4')
        end

        
        
        %%X(204, 24), y(6,24)
        function showCorrelationEachBand(X, y)
            lenY = size(y, 1);
            
            figure
%             tiledlayout(2, 3);
            hold on
            
            for i = 1 : lenY
               yy = y(i, :);
               
               lenX = size(X, 1);
               corss = zeros(lenX, 1);
               
               for j = 1 : lenX
                  xx = X(j, :);
                  tmp = corrcoef(xx, yy);
                  corss(j) = tmp(1,2);
               end
%                nexttile
               plot(corss);
            end            
            hold off
            legend('As', 'Cd', 'Cr', 'Cu', 'Ni', 'Pb')
            
% %             A = corrcoef(X');
% %             figure
% %             h = heatmap(A);

        end
        
        
        function plotPredict(y, yfit)            
            lenY = size(y, 2);
            
            figure
            tiledlayout(2, 3);
%             tiledlayout('flow');
            
            for i = 1 : lenY
               yy = y(:, i);
               yhat = yfit(:, i);
               
               nexttile
               hold on
               plot(yy, yhat, 'o');
               ymin = min(min(yy, [], 'all'), min(yhat, [], 'all'));
               ymax = max(max(yy, [], 'all'), max(yhat, [], 'all'));
               
               
               dd = (ymax - ymin) / 23;
               kkk = ymin * 0.9: dd: ymax*1.1;
               plot(kkk, kkk, '--k');         
               
% %                for j = 1 : length(yy)                   
% %                   text(yy(j)*0.97, yhat(j)*1.05, num2str(j)) 
% %                end
               hold off
            end            
        end
        %%画出每个样本的两个采样曲线
        function plotPair(XX)
            colors = [[1, 0, 0]; [0, 1, 0]; [0, 0, 1]; [0, 0, 0]; 
                      [1, 0, 1]; [0, 1, 1]; [0.5, 0.5, 0.5]; [0.7, 0.2, 0.2];
                      [0.3, 0.7, 0.1]; [0.3, 0.2, 0.7]; [0.7, 0.7, 0.1]; [0.7, 0.2, 0.7];
                      ];
                  
            lenX = size(XX, 2)/2;
            figure;
            hold on
            for i = 1 : lenX
                plot(XX(:, i), 'Color', colors(i, :));                
                plot(XX(:, i+lenX), 'Color', colors(i, :));
            end
            legend('1', '1', '2', '2', '3', '3', '4', '4', '5', '5', '6', '6')
            hold off
        end
        
        
        function plotTest()
            x=-3:.25:5;
            y=-0.3*x+3.5.*x.^2-x.^3+20*rand([1,length(x)]);
            [p,S]=polyfit(x,y,3);
            % 计算以p为系数的多项式在 x 中各点处的拟合值。将误差估计结构体指定为第三个输入，
            % 以便polyval 计算标准误差的估计值。标准误差估计值在 delta 中返回。
            
            [y_fit,delta]=polyval(p,x,S);
            % 绘制原始数据、线性拟合和 95% 预测区间 y±2Δ。
            uy=y_fit+2*delta;
            dy=y_fit-2*delta;
            % 绘制原始数据
            plot(x,y,'rx','LineWidth',1.2)
            hold on
            % 绘制拟合曲线
            plot(x,y_fit,'Color',[82,124,179]./255,'LineWidth',1.5)
            % 绘制置信区间
            plot([x',x'],[uy',dy'],'Color',[82,124,179]./255,'LineWidth',1.2,'LineStyle','--')
            fill([x,x(end:-1:1)],[uy,dy(end:-1:1)],[82,124,179]./255,'EdgeColor','none','FaceAlpha',.2)
            
            title('Linear Fit of Data with 95% Prediction Interval')
            legend('Data','Linear Fit','95% Prediction Interval')
            
        end
        function plotRegion(x, y, std)
            figure
            hold on
            plot(x, y, 'r');
            uy = y + std;
            dy = y - std;
            
            plot(x, uy, 'b--')
            plot(x, dy, 'b--')
            fill([x;x(end:-1:1)], [uy;dy(end:-1:1)],[82,124,179]./255,'EdgeColor','none','FaceAlpha',.2)
            
            hold off
            
        end
        
        %%% 画出12条光谱曲线
        function plot12Curve(xx, sTitle)
            figure
            hold on
            for i = 1 : length(xx)
                x = xx{i};
                plot(x(:, 1), x(:, 4))
            end
            legend()
            title(sTitle)
            hold off
        end
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        function x_data =  extractMatrix(xx, comp)
            sz1 = size(xx{1},1);
            
            x_data = zeros(sz1, length(xx));
            for i = 1 : length(xx)
                x = xx{i};
                
                x_data(:, i) = x(:, comp);
            end
        end
        %%添加X轴标签
        function AddXTickLabels()
%             xticklabels({'ph值', '有机质', '有效磷', '速效钾', '全氮', '砷', '镉', '铬', '铜', '镍', '铅'})
            xticklabels({'砷', '镉', '铬', '铜', '镍', '铅'})
        end
        
        function AddLegends()
%             legend('ph值', '有机质', '有效磷', '速效钾', '全氮', '砷', '镉', '铬', '铜', '镍', '铅')
            legend('砷', '镉', '铬', '铜', '镍', '铅')
        end
        
        %         function SaveStatsToFiles(fileSave, rmse, R2, R2total, rsd, rpd)
        %             data = [rmse; R2; rsd; rpd];
        %             save(fileSave, '-ascii', '-tabs', 'data');
        %             save(fileSave, '-ascii', '-tabs', '-append', 'R2total');
        %         end
        
        function SaveStatsToFiles(fileSave, arg)
            % % % 代码OK，但要增加MAPE，不需要R2total
            % % %             data = [arg.rmse; arg.R2; arg.rsd; arg.rpd];
            % % %             save(fileSave, '-ascii', '-tabs', 'data');
            % % %             R2total = arg.R2total;
            % % %             save(fileSave, '-ascii', '-tabs', '-append', 'R2total');
            
            data = [arg.rmse; arg.R2; arg.rsd; arg.rpd; arg.mape];
            save(fileSave, '-ascii', '-tabs', 'data');
        end
        
        function SaveValuesToFile(fileSave, arg)
           data =  [arg.y, arg.yfit];
           save(fileSave, '-ascii', '-tabs', 'data');
        end
        
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        % % % % % % % % % % % % % % % % % % % % % % % % % % % %         %%
        function ret = CalcAllErrors(y, yfit, bShow, cc)
            if nargin < 4
                ret = Helper.CalcAllErrorsPatch(y, yfit, bShow);
                ret.y = y;
                ret.yfit = yfit;
                return
            end
            %%只针对 24折进行处理， 不是24折的，暂不处理
            if cc.NumTestSets ~= 24
               disp('error in 24-fold') 
            end
  
        end
        
        %% 想改但有没改成 2024-02-02
        function ret = CalcAllErrorsPatch(y, yfit, bShow)
             %% y =(24, 6)
             nn = size(y, 1);
             meanY = mean(y);
             rmse = rms(y - yfit);
             rmse_std = std(y - yfit);
             rsd = rmse ./ meanY;%%标准偏差系数 relative standard deviation
             %% rsd 应该给错了公式 2024-02-01， 应该是标准差/平均值
             
             
             %            rmse = rmse * sqrt(nn / (nn -1));
             
             %             R2total = norm(yfit - meanY)^2 / norm(y - meanY)^2;
             %             R2  = vecnorm(yfit - meanY).^2 ./ vecnorm(y - meanY).^2;
             R2total = 1 - norm(y - yfit)^2 / norm(y - meanY)^2;
             R2  = 1 - vecnorm(y - yfit).^2 ./ vecnorm(y - meanY).^2;
             
             ss = std(y); %% 分母为N-1，
             %             ss = std(y,1); %% 2024-01-30 改为分母N，
             rpd = ss ./ rmse; %%Relative Percentage Difference  residual prediction deviation
             
             
             err = abs(y - yfit);
             mape = mean(err ./ y);
             mape_std = std(err ./ y);
             
             if bShow
                 figure
                 hold on
                 xx = 1: size(y, 2);
                 plot(xx, rsd, '-rs')
                 plot(xx, R2, '-bo')
                 %                plot(rpd, '-kd')
                 hold off
                 %            legend('ph值', '有机质', '有效磷', '速效钾', '全氮', '砷', '镉', '铬', '铜', '镍', '铅')
                 xticks(xx)
                 Helper.AddXTickLabels()
                 title(R2total)
                 legend('相对标准偏差', 'R2')
             end
             ret.rmse = rmse;
             ret.R2 = R2;
             ret.R2total = R2total;
             ret.rsd = rsd;
             ret.rpd = rpd;
             ret.mape = mape;
             ret.mape_std = mape_std;
             ret.rmse_std = rmse_std;
        end
       
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        function PlotStemErrors(y, y_hat)
            residuals = y - y_hat;
            figure
            
            stem(residuals)
            xlabel('Observation');
            ylabel('Residual');
            legend
            
        end
        
        function PlotPLSRComps(X, y, nn)
            [XL,YL,XS,YS, beta, PCTVAR] = plsregress(X, y, nn);
            Helper.PlotPLSRvar(PCTVAR, true);
        end
        %%画出Y的累计方差
        function PlotPLSRvar(PCTVAR, bNewFig)
            nn = size(PCTVAR, 2);
            if bNewFig == true
                figure
            end
            hold on
            %             plot(1:nn,cumsum(100*PCTVAR(2,:)),'-bo');
            plot(1:nn,cumsum(100*PCTVAR(2,:)));
            xlabel('Number of PLS components');
            ylabel('Percent Variance Explained in y');
            
            %             plot(1:nn,cumsum(100*PCTVAR(1,:)),'-rs');            
            hold off
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%设置cv-12，使每对取样一致 2023-08-08
        function cc = GetMyCVParition()
            cc = mycvpartition(24, 'KFold', 12);%%12折
%             for i = 1 : 12
%                 trai2 = true(24, 1);
%                 
%                 trai2(i) = 0;
%                 trai2(12 + i) = 0;
%                 
%                 test2 = ~trai2;
%                 cc.training(i) = trai2;
%                 cc.test(i) = test2;
%             end
        end
       
        function ret = isModelling(step)
            step4 = 8;
            step5 = 16;     step6 = 32;    step7 = 128;     step8 = 256;

           ret =  bitand(step, step4 | step5 | step6 | step7 | step8);
        end
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %         
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
        function ret = DoAllMetalsModelFeaturePreProcess(x_data, y_data, arg)
               
            methods = arg.methods;
            ret = cell(size(methods, 1), 1);
            
            for metal = 1 : size(methods, 1)%%元素个数
                
                preType = methods(metal, 1);
                modelType = methods(metal, 2);
                
                if preType == 0 && modelType == 0
                    disp('preType == 0 & modelType == 0')
                    continue; 
                end
                ret{metal} = Helper.DoOneMetalWithFea(x_data, y_data, metal, preType, modelType, arg);              
            end
        end

        function ret = DoOneMetalWithFea(X, y, metal, preType, modelType, arg)
            %%X，24*204， y，24*6       
            bShow = arg.bShow;            
        
            x_data = PreHelper.DoOnePreProcessing(X, preType, false);          
            
            arg.metal = metal;
            arg.wave = 10;%%10个波段选择
            arg.bEach = false;            
            
            ret = cell(arg.feaTypes, 1);            
            
            % if modelType ~= 4 || metal ~= 3
            %     return;
            % end
            feaTypes = arg.feaTypes;

            for feaType = 1 : feaTypes %2024-07-27 %%专门针对 3-5的特征方法进行
                disp(['----fea：', num2str(feaType), '----model:', num2str(modelType), '----metal:', num2str(metal)]);
                arg.feaType = feaType;
                ret{feaType} = MyModel.DoOneModel(x_data', y', modelType, bShow, arg);
            end
        end
        %%2024-05-01 设置的每一个元素与指定settting，返回opt
        function ret = DoOneMetalWithOneSetting(X, y, metal, preType, feaType, arg)
            bShow = arg.bShow;
            
            x_data = PreHelper.DoOnePreProcessing(X, preType, false);
            arg.metal = metal;
            arg.wave = arg.wave;
            arg.bEach = false;
            arg.ncomp = 10;%%%%%%%%有问题要改
            
            disp(['----fea：', num2str(feaType), '----preType:', num2str(preType), '----metal:', num2str(metal)]);
            
            kopt = FeaSelect.GetOptimalFeature(x_data', y', arg.wave, feaType, metal, bShow);
            ret = sort(kopt);
% %             ret = kopt{metal};
% %             ret = sort(ret);
        end
        
        function data = SetResultForOneMetal(ret)
            len = length(ret);%%=6 金属个数
            data = cell(len, 1);
            for i = 1 : len
                metal = ret{i};
                if isempty(metal) == true
                    disp(['metal ', num2str(i), ' is empty'])
                    continue;
                end
                
                sz = length(metal);%%特征方法个数
                
                data{i}.rmse = zeros(1, sz);
                data{i}.R2 = zeros(1, sz);
                data{i}.rsd = zeros(1, sz);
                data{i}.rpd = zeros(1, sz);
                data{i}.mape = zeros(1, sz);
                data{i}.mape_std = zeros(1, sz);
                data{i}.rmse_std = zeros(1, sz);
                for j = 1 : sz
                    data{i}.rmse(j) = metal{j}.rmse(i);
                    data{i}.R2(j) = metal{j}.R2(i);
                    data{i}.rsd(j) = metal{j}.rsd(i);
                    data{i}.rpd(j) = metal{j}.rpd(i);
                    data{i}.mape(j) = metal{j}.mape(i);
                    data{i}.mape_std(j) = metal{j}.mape_std(i);
                    data{i}.rmse_std(j) = metal{j}.rmse_std(i);
                end
            end
        end
        % % % % % % % % % % % % % % % % % % % % % % % % % %
        function ret = DoOneModelWithAllSetting(X, y, modelType, arg)
            % % 2024-08-05 一个模型，一个元素            
            
            arg.bEach = false;
            
            ret = cell(arg.preTypes,1);
            
            bWait = isfield(arg, 'wait');%%2024-02-10增加的特征选择！
            for preType = 1 : arg.preTypes
                
                x_data = PreHelper.DoOnePreProcessing(X, preType, false);                

                disp(['----preType：', num2str(preType), '----modelType:', num2str(modelType)]);
                ret{preType} = MyModel.DoOneModel(x_data', y', modelType, arg.bShow, arg);   
                
                if bWait == true               
                    pause(arg.wait) 
                end
            end            
        end
        
        
        
        % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %
        %% 在平均光谱线上 画出特征波段 2024-07-26
        function DrawSelectionOnMeanSpectral(x_data, fsss, offStart, offEnd)
            data = mean(x_data);

            lenX = length(data) + offStart + offEnd;
            for i = 1:length(fsss)
                xx = fsss{i}+offStart;
                yy = data(xx);
                
                figure
                plot(1:lenX, data);
                hold on
                plot(xx, yy, 'o');
                hold off
            end
            xlim([1, lenX]);
            
        end

        %%% x_data  24*206
        function fig = DrawSelectionOnMeanSpectral2(x_data, wave, offStart, offEnd, bandX)
            data = mean(x_data);
            
            lenX = length(data) + offStart + offEnd;
            
            pos = wave + offStart
            xx = bandX(pos);
            yy = data(pos);
            
            fig = figure
            % %            plot(1:lenX, data);
            plot(bandX, data);
            
            hold on
            plot(xx, yy, 'o', 'MarkerSize', 8, 'LineWidth',2);
            hold off
            % %             xlim([1, lenX]);
            xlim([min(bandX)-3, max(bandX)+3]);
            ylim([0.38, 0.44])
            xlabel('Wavelength (nm)')
            ylabel('Reflentance');
        end
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % %         
        %% 用作合并两个输出文件 2024-07-29
        function ret = MergeRetFiles(file1, file2, fileSave, nCut)
            ret1 = load(file1, 'ret');
            ret2 = load(file2, 'ret');
            ret1 = ret1.ret;
            ret2 = ret2.ret;
            
            if length(ret1) ~= length(ret2)
                disp('error in length')
                return
            end            
            ret = cell(length(ret1), 1);
            for i = 1 : length(ret1)
                data1 = ret1{i};
                data2 = ret2{i};
                ret{i} = [data1(1:nCut-1); data2(nCut : end)];
            end
            
            save(fileSave, 'ret');
        end
    end
end
    
   